Los datos obtenidos de Kaggle tienen la información atmosfĆ©rica de varios aƱos en Australia, datos en los que viene incluida desde la velocidad del viento hasta la temperatura, en total 23 variables recogidas durante varios aƱos en distintas ciudades de Australia con las que se obtienen un dataset de 140.000 lĆneas aproximadamente.
Nuestro objetivo serĆ” predecir la lluvia del dĆa siguiente con los datos metereológicos del dĆa.
Con este dataset tan grande lo primero que nos planteamos fue centrarnos en dos cosas:
Utilizar una zona concreta de Australia, sacada de la variable localización, de la cual elegimos 4 ciudades situadas en la costa sureste de Australia.
Utilizar la variable temporal de alguna forma, ya que considerabamos que tenĆa importancia pero no podĆamos usar cada dĆa del aƱo como un dato diferente, por lo que decidimos obtener a partir de la fecha la estación del aƱo en la que estaba cada lĆnea.
## Dimensiones dataset train: 9824 25
## Dimensiones dataset test: 1228 25
## Dimensiones dataset validación: 1228 25
Antes que nada, visualizamos información bÔsica de las ciudades elegidas y estaciones y cómo estÔn relacionadas con la lluvia.
Realizamos un conteo del nĆŗmero de dĆas que han llovido o no en cada una.
Si filtramos por los dĆas en los que ha llovido hoy (RainToday = 1), vemos las veces que ha llovido el dĆa siguiente o no.
Vemos que en Adelaide, Canberra y Melbourne no suele ser fuecuente que llueva el dĆa siguiente si ha llovido hoy. En Sydney en cambio, es mĆ”s fequente que llueva si ha llovido el dĆa anterior.
Representamos ahora los dĆas que han llovido en función de las ciudades y las estaciones del aƱo.
En Adelaide y Melbourne las estaciones sĆ influyen mĆ”s en la frecuencia de dĆas que llueven, mientras que en Sydney y Canberra suele ser mĆ”s homogĆ©nero.
Resumen de las variables.
## Date Season Location
## Min. :2007-11-01 Length:9824 Length:9824
## 1st Qu.:2010-05-14 Class :character Class :character
## Median :2012-08-31 Mode :character Mode :character
## Mean :2012-10-06
## 3rd Qu.:2015-01-23
## Max. :2017-06-25
##
## MinTemp MaxTemp Rainfall Evaporation
## Min. :-8.00 Min. : 4.10 Min. : 0.000 Min. : 0.000
## 1st Qu.: 7.90 1st Qu.:17.10 1st Qu.: 0.000 1st Qu.: 2.400
## Median :11.60 Median :21.40 Median : 0.000 Median : 4.200
## Mean :11.44 Mean :22.01 Mean : 2.165 Mean : 5.018
## 3rd Qu.:15.50 3rd Qu.:26.10 3rd Qu.: 0.600 3rd Qu.: 6.800
## Max. :33.90 Max. :45.80 Max. :119.400 Max. :43.400
## NA's :9 NA's :8 NA's :181 NA's :2493
## Sunshine WindGustDir WindGustSpeed WindDir9am
## Min. : 0.00 Length:9824 Min. : 11.00 Length:9824
## 1st Qu.: 4.10 Class :character 1st Qu.: 31.00 Class :character
## Median : 7.90 Mode :character Median : 39.00 Mode :character
## Mean : 7.15 Mean : 40.64
## 3rd Qu.:10.20 3rd Qu.: 48.00
## Max. :13.90 Max. :106.00
## NA's :2667 NA's :1104
## WindDir3pm WindSpeed9am WindSpeed3pm Humidity9am
## Length:9824 Min. : 0.00 Min. : 0.00 Min. : 5.00
## Class :character 1st Qu.: 7.00 1st Qu.:13.00 1st Qu.: 56.00
## Mode :character Median :11.00 Median :19.00 Median : 68.00
## Mean :13.27 Mean :18.82 Mean : 67.22
## 3rd Qu.:19.00 3rd Qu.:24.00 3rd Qu.: 80.00
## Max. :63.00 Max. :65.00 Max. :100.00
## NA's :207 NA's :199 NA's :71
## Humidity3pm Pressure9am Pressure3pm Cloud9am
## Min. : 4.00 Min. : 982.3 Min. : 985.5 Min. :0.000
## 1st Qu.:37.00 1st Qu.:1013.6 1st Qu.:1011.4 1st Qu.:2.000
## Median :48.00 Median :1018.6 Median :1016.3 Median :6.000
## Mean :48.95 Mean :1018.5 Mean :1016.2 Mean :4.694
## 3rd Qu.:60.00 3rd Qu.:1023.5 3rd Qu.:1021.1 3rd Qu.:7.000
## Max. :99.00 Max. :1040.2 Max. :1037.8 Max. :9.000
## NA's :30 NA's :201 NA's :196 NA's :4111
## Cloud3pm Temp9am Temp3pm RainToday
## Min. :0.000 Min. :-1.30 Min. : 3.70 Length:9824
## 1st Qu.:2.000 1st Qu.:11.70 1st Qu.:15.90 Class :character
## Median :5.000 Median :15.40 Median :20.00 Mode :character
## Mean :4.696 Mean :15.55 Mean :20.56
## 3rd Qu.:7.000 3rd Qu.:19.30 3rd Qu.:24.50
## Max. :8.000 Max. :38.60 Max. :44.70
## NA's :4298 NA's :24 NA's :17
## RISK_MM RainTomorrow
## Min. : 0.000 Length:9824
## 1st Qu.: 0.000 Class :character
## Median : 0.000 Mode :character
## Mean : 2.164
## 3rd Qu.: 0.800
## Max. :119.400
##
Analizamos las variables individuales por separado con distintos grƔficos.
En este dataset hay muchos pares de variables que estĆ”n fuertemente relacionadas, por ejemplo la temperatura mĆ”xima y mĆnima de un dĆa, o la presión a las 9 de la maƱana y la presión a las 3 de la tarde. Por ello, en el anĆ”lisis individual de variables se estudiarĆ”n a la vez por una mejor comprensión.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
La temperatura mĆnima se podrĆa asimilar a una distribución normal. La temperatura mĆ”xima tiene una cola a su derecha en la que aparecen mĆŗtliples valores atĆpicos.
Las variables Temp9am y Temp3pm son muy parecidas a las temperaturas mĆnimas y mĆ”ximas respectivamente.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Estas variables tienen un comportamiento parecido a la temperatura mĆnima y mĆ”xima respectivamente.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Parecen tener una distribución normal ambas variables.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Se trata de una variable categórica que indica la fracción del cielo que estĆ” cubierta por nubes a las 9am. EstĆ” medida en āoktasā, que son unidades de octavos. Una medida cero indica que el cielo estĆ” completamente despejado y un ocho indica que estĆ” completamente cubierto.
La variable a las diferentes horas se distribuye de forma muy parecida.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
La velocidad del viente suele tomar valores mayores por la tarde, a las 3pm comparado con la maƱana.
Las direcciones en las que el viento se mueve son muy diferentes a las 9 y a las 3. Predominan la dirección Norte y la Oeste por la mañana, mientras que por la tarde suele ser mÔs homogéneo.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Se trata de la dirección y la velocidad del viento mÔs fuerte en las últimas 24 horas. Predomina el viento en la dirección norte con una media de veolcidad de 40km/h.
Analizamos cuatro variables que no estƔn, a priori, relacionadas por pares.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Es la cantidad de lluvia que se ha registrado en el dĆa en mm. Como la mayorĆa de dĆas no llueve, aparecen muchos valores nulos en el histograma que dificulta la normalización de la variable.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Se trata de la evaporación en mm en las Ćŗltimas 24h a las 9am. Contiene un gran nĆŗmero de valroes atĆpicos.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Sorprende el nĆŗmero de dĆas en el que las horas de sol es 0, es decir, estĆ” nublado todo el dĆa o con rayos de sol muy dĆ©biles.
Analizamos la relación unas variables con otras.
Gracias a este grÔfico demostramos que los pares de variables citados sà estÔn fuertemente relacionados. Por ejemplo, la presión a las 9 de la mañana con la presión a las 3 de la tarde: si una aumenta, la otra también. Destacar que también hay correlaciones inversas: cuando aumenta la variable Sunshine, disminuye la Cloud9am.
Analizamos en general la relación entre las variables con las estaciones, la variable de salida (RainTomorrow) y las ciudades. Para ello seleccionamos una de las variables de los pares (ya que contienen información similar) y el resto. La variable Rainfall no la mostramos pues su distribución es difĆcil de visualizar. Se analizarĆ”n posteriormente con las transformaciones.
Relaciones por estaciones.
La Temperatura (MaxTemp, Temp9am), la presión (Pressure3pm) y evaporación muestran un claro comportamiento diferente según la estación.
Relaciones por ciudades.
La temperatura es la variable dónde se puede observar mÔs claramente que tiene un comportamiento diferente para cada ciudad.
Relaciones por RainTomorrow, si llueve o no.
Las variables interesantes a analizar con mÔs detalle son: la temperatura, la humedad, la presión y los rayos de sol.
A partir de estas relaciones, indagamos con mƔs detalle las relaciones que parecen interesantes.
Comprobamos que, efectivamente, las cuatro variables de temperatura son muy parecidas, como puede observarse en sus distribuciones:
Al estar relacionadas y tener un comportamiento similar posteriormente se estudiarƔ introdudir al modelo interacciones entre Ʃstas.
Veamos cómo se comporta una de ellas según las estaciones y ciudades:
La Temperatura claramente estĆ” condicionada por las estaciones y por las ciudades. No hay mucha diferencia entre otoƱo y primavera. Dada esta información se podrĆa plantear realizar un modelo para cada estación o para cada ciudad.
Comprobamos que los pares de temperaturas estƔn fuertemente relacionados y que por lo tanto serƔ interesante estudiar interacciones entre ellas.
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Analizamos otras variables mÔs en detalle. Veamos cómo se comporta la presión según las estaciones y ciudades:
La presión se diferencia bastante con cada estación y también se aprecia que estÔ bastante relacionada con si llueve o no: la presión suele ser menor cuando llueve.
Comprobamos que las variables de presión estÔn fuertemente relacionadas entre ellas.
Veamos cómo se comporta la humedad según las estaciones y ciudades:
Vemos que la humedad contiene valores diferentes segĆŗn si llueve o no y se diferencia bĆ”sicamente entre otoƱo e invierno. Efectiavmente la humedad es mayor los dĆas que llueve.
Comprobamos cómo de relacionadas estÔn las variables de humedad entre ellas.
Veamos cómo se comporta la evaporación según las estaciones y ciudades:
Se observan comportamientos diferentes de la evaporación sobretodo en función de las estaciones.
Veamos cómo se comportan las horas de sol según las estaciones y ciudades:
Claramente la variable sunshine se diferencia mÔs en función de la variable de salida RainTomorrow.
Veamos cómo se comporta la velocidad del viento según las estaciones y ciudades:
La velocidad del viento suele ser mayor los dĆas que llueve.
Veamos cómo se comportan las nubes del viento según las estaciones y ciudades:
Las variables de nubes a las diferentes horas se comportan de manera muy parecida.
Veamos cómo se comportan la velocidad del viento y su dirección según las estaciones y ciudades:
No parece haber mucha diferencia entre la dirección del viento y la velocidad en función de si llueve o no.
AƱadir estos grƔficos (Rainfall) despuƩs de analizar las transformaciones.
Se puede observar en la grafica siguiente que ninguna variable individual es mejor que otra para diferenciar si llueve o no llueve, estudiaremes MinTemp, MaxTemp, Temp9am, Temp3pm, Pressure9am, Pressure3pm, Humidity9am, Humidity3pm.
Ahora crearemos unas nuevas varibles combinandolas entre si. Las usaremos en el modelo mas adelante para ajustarlo.
Vamos a realizar un estudio de los datos faltantes para despuƩs imputarlos utilizando las tecnicas junto al resto de parametros.
Para analizar los posibles datos faltantes y las relaciones utilizamos dos bibliotecas, visdat y naniar. Con la base de datos de train vemos como se reparten los datos faltantes
## Date Season Location MinTemp MaxTemp
## 0 0 0 9 8
## Rainfall Evaporation Sunshine WindGustDir WindGustSpeed
## 181 2493 2667 1105 1104
## WindDir9am WindDir3pm WindSpeed9am WindSpeed3pm Humidity9am
## 709 229 207 199 71
## Humidity3pm Pressure9am Pressure3pm Cloud9am Cloud3pm
## 30 201 196 4111 4298
## Temp9am Temp3pm RainToday RISK_MM RainTomorrow
## 24 17 181 0 0
## AvgTemp Temp Pressure Humidity
## 4 2 178 9
Con estos Datos podemos deducir que las variables con mas datos faltantes serian (de mas a menos datos faltantes): - Cloud3pm - Cloud9am - Sunshine - Evaporation - WindGustDir - WindGustSpeed - WindDir9am
Con estos datos vamos a estudiar cuales son las variables con mas relaciones de los valores faltantes
Con esta grafica podemos observar que los datos faltantes mas relaccionados son sobretodo Cloud3pm y Cloud9am, esto puede deberse a que estas variables son categoricas y no pueden medirse directamente, ya que necesitan de la evaluacion del cielo en ese dia y esto puede que este automatizado.
A estas variables se le suman las variables Evaporation y Sunshine, que tambien estan relaccionadas tanto entre si como con las mencionadas anteriormente.
Por ultimo recalcar que la Variable WindGustDir tiene muchos datos faltantes en los que no tiene relacion con las antes mencionadas. Esto puede deberse a que miden cosas completamente diferentes y no necesariamente dependen de las otras.
Con todos estos datos y el EDA con el que hemos empezado podemos comenzar a imputar datos faltantes:
Las primeras variables que vamos a evaluar para completar sus datos faltantes son aquellas variables numericas que se miden en dos momentos del dia, esto nos puede ayudar para completar una con la otra ya que en gran medida estas variables son bastannte dependientes y podemos sacar una buena correlacion para completar una con la otra.
Para esto utilizaremos el metodo de imputacion KNN, con estos datos consideramos que se pueden completar perfectamente estos valores faltantes.
Variables Temp9am & Temp3pm:
#Variables Temp9am & Temp3pm: método: imputación simple: KNN
par(mfrow=c(1,1))
train %>% select(Temp9am, Temp3pm) %>% marginplot()Viendo esta grafica podemos decir que tiene sentido utilizar la imputacion mediante KNN.
#Posible imputación simple mediante un algoritmo de clustering k-NN (con la función VIM:kNN)
#Parece que tiene sentido:
train %>% select(Temp9am, Temp3pm) %>% VIM::kNN() %>% marginplot(., delimiter="_imp")train_imputed <- kNN(train, variable=c("Temp9am","Temp3pm"))
#test_imputed <- kNN(test, variable=c("Temp9am","Temp3pm"))Comprobamos si nuestra muestra inicial modifica su densidad para asegurarnos de que esta imputacion funciona y no hemos generado ningun valor que no coincida.
#Comprobaciones:
par(mfrow=c(1,2))
plot(density(train$Temp9am, na.rm = T), col=2, main="Temp9am")
lines(density(train_imputed$Temp9am), col=3)
plot(density(train$Temp3pm, na.rm = T), col=2, main="Temp3pm")
lines(density(train_imputed$Temp3pm), col=3) Como podemos comprobar, apenas se aprecia la linea roja que esta debajo de la densidad inicial por lo que podemos decir que no cambia
Variables MinTemp & MaxTemp:
#Variables MinTemp & MaxTemp: método: imputación simple: KNN
par(mfrow=c(1,1))
train %>% select(MinTemp, MaxTemp) %>% marginplot()#Parece que tiene sentido:
train %>% select(MinTemp, MaxTemp) %>% VIM::kNN() %>% marginplot(., delimiter="_imp")train_imputed <- kNN(train, variable=c("MaxTemp","MinTemp"))
#test_imputed <- kNN(test, variable=c("MaxTemp","MinTemp"))#Comprobaciones:
par(mfrow=c(1,2))
plot(density(train$MaxTemp,na.rm = T),col=2,main="MaxTemp")
lines(density(train_imputed$MaxTemp),col=3)
plot(density(train$MinTemp,na.rm = T),col=2,main="MinTemp")
lines(density(train_imputed$MinTemp),col=3) En este caso la imputacion de datos faltantes tambien funciona
Variables Pressure9am & Pressure3pm:
#Variables Pressure3pm & Pressure9am: método: imputación simple: KNN
par(mfrow=c(1,1))
train %>% select(Pressure3pm, Pressure9am) %>% marginplot()#Parece que tiene sentido:
train %>% select(Pressure3pm, Pressure9am) %>% VIM::kNN() %>% marginplot(., delimiter="_imp")train_imputed <- kNN(train_imputed, variable=c("Pressure3pm","Pressure9am"))
#test_imputed <- kNN(test_imputed, variable=c("Pressure3pm","Pressure9am"))#Comprobaciones:
par(mfrow=c(1,2))
plot(density(train$Pressure3pm,na.rm = T),col=2,main="Pressure3pm")
lines(density(train_imputed$Pressure3pm),col=3)
plot(density(train$Pressure9am,na.rm = T),col=2,main="Pressure9am")
lines(density(train_imputed$Pressure9am),col=3) Podemos utilizar esta imputacion en la presion.
Variables Humidity9am & Humidity3pm:
#Variables Humidity9am & Humidity3pm: método: imputación simple: KNN
par(mfrow=c(1,1))
train %>% select(Humidity3pm, Humidity9am) %>% marginplot()#Parece que tiene sentido:
train %>% select(Humidity3pm, Humidity9am) %>% VIM::kNN() %>% marginplot(., delimiter="_imp")train_imputed <- kNN(train_imputed, variable=c("Humidity3pm","Humidity9am"))
#test_imputed <- kNN(test_imputed, variable=c("Humidity3pm","Humidity9am"))#Comprobaciones:
par(mfrow=c(1,2))
plot(density(train$Humidity3pm,na.rm = T),col=2,main="Humidity3pm")
lines(density(train_imputed$Humidity3pm),col=3)
plot(density(train$Humidity9am,na.rm = T),col=2,main="Humidity9am")
lines(density(train_imputed$Humidity9am),col=3) Podemos utilizar esta imputacion en la humedad.
Variables Rainfall & Raintoday: Podriamos pensar que estas dos variables van de la mano, ya que si la premisa se cumple, cuando la variable Rainfall es 0 es que ese mismo dia no ha llovido, por ello haremos un estudio de los datos faltantes.
Esto quiere decir que los datos faltantes son exactamente los mismos, por lo que vamos a proceder a completar los datos faltantes de la variable rainfall y con ellos completaremos los datos de Raintoday.
Otras variables que puede tener sentido relacionar es la cantidad de lluvia con la humedad.
Variables Humidity9am & Rainfall:
#Variables Humidity9am & Rainfall: método: imputación simple: KNN
par(mfrow=c(1,1))
train %>% select(Humidity9am, Rainfall) %>% marginplot()#Parece que tiene sentido:
train %>% select(Humidity9am, Rainfall) %>% VIM::kNN() %>% marginplot(., delimiter="_imp")train_imputed <- kNN(train_imputed, variable=c("Humidity9am","Rainfall"))
# test_imputed <- kNN(test_imputed, variable=c("Humidity9am","Rainfall"))train_imputed[c("Temp9am_imp","Temp3pm_imp" ,"MaxTemp_imp","MinTemp_imp", "Humidity9am_imp","Rainfall_imp", "Humidity3pm_imp", "Pressure3pm_imp","Pressure9am_imp")] <- list(NULL)
imputacion_knn_test<-function(train_imputed, test){
train_imputed["dataset"] <- 'train'
test["dataset"] <- 'test'
dataset_knn = rbind(train_imputed, test)
dataset_knn_imputed <- kNN(dataset_knn, variable=c("Temp9am","Temp3pm"))
dataset_knn_imputed <- kNN(dataset_knn, variable=c("MaxTemp","MinTemp"))
dataset_knn_imputed <- kNN(dataset_knn, variable=c("Humidity3pm","Humidity9am"))
dataset_knn_imputed <- kNN(dataset_knn, variable=c("Pressure3pm","Pressure9am"))
dataset_knn_imputed <- kNN(dataset_knn, variable=c("Humidity9am","Rainfall"))
train_imputed = subset(dataset_knn_imputed, dataset == 'train')
test_imputed = subset(dataset_knn_imputed, dataset == 'test')
test_imputed[c("Temp9am_imp","Temp3pm_imp" ,"MaxTemp_imp","MinTemp_imp", "Humidity9am_imp","Rainfall_imp", "Humidity3pm_imp", "Pressure3pm_imp","Pressure9am_imp","dataset")] <- list(NULL)
return (test_imputed) }
test_imputed = imputacion_knn_test(train_imputed, test)#Comprobaciones:
par(mfrow=c(1,2))
plot(density(train$Humidity9am,na.rm = T),col=2,main="Humidity9am")
lines(density(train_imputed$Humidity9am),col=3)
plot(density(train$Rainfall,na.rm = T),col=2,main="Rainfall")
lines(density(train_imputed$Rainfall),col=3)## Date Season Location MinTemp MaxTemp
## 0 0 0 0 0
## Rainfall Evaporation Sunshine WindGustDir WindGustSpeed
## 0 2493 2667 1105 1104
## WindDir9am WindDir3pm WindSpeed9am WindSpeed3pm Humidity9am
## 709 229 207 199 0
## Humidity3pm Pressure9am Pressure3pm Cloud9am Cloud3pm
## 0 0 0 4111 4298
## Temp9am Temp3pm RainToday RISK_MM RainTomorrow
## 24 17 181 0 0
## AvgTemp Temp Pressure Humidity
## 4 2 178 9
Ahora que ya hemos logrado imputar la variable de Rainfall la utilizaremos para conseguir los datos faltantes de raintoday:
train_imputed$RainToday[which(is.na(train_imputed$RainToday))] <- ifelse(train_imputed$Rainfall == 0, "No", "Yes")
test_imputed$RainToday[which(is.na(test_imputed$RainToday))] <- ifelse(test_imputed$Rainfall == 0, "No", "Yes")
train_imputed$RainToday <- ifelse(train_imputed$RainToday == "Yes", 1, 0)
train_imputed$RainToday = as.numeric(train_imputed$RainToday)
train_imputed$RainTomorrow <- ifelse(train_imputed$RainTomorrow == "Yes", 1, 0)
train_imputed$RainTomorrow = as.numeric(train_imputed$RainTomorrow )
test_imputed$RainToday <- ifelse(test_imputed$RainToday == "Yes", 1, 0)
test_imputed$RainToday = as.numeric(test_imputed$RainToday)
test_imputed$RainTomorrow <- ifelse(test_imputed$RainTomorrow == "Yes", 1, 0)
test_imputed$RainTomorrow = as.numeric(test_imputed$RainTomorrow )Variables Cloud3pm & Cloud9am: Para estas variable sutilizamos la libreria mice con el metodo de la imputacion multiple y comprobamos que no varia la densidad de los valores
par(mfrow=c(1,1))
train_imputed0 <- mice(train_imputed[,c('Cloud3pm','Cloud9am')], seed=2018, print = F, m = 30)
train_imputed1 <- mice::complete(train_imputed0)
test_imputed0 <- mice(test_imputed[,c('Cloud3pm','Cloud9am')], seed=2018, print = F, m = 30)
test_imputed1 <- mice::complete(test_imputed0)
xyplot(train_imputed0, Cloud3pm ~Cloud9am)par(mfrow=c(1,2))
plot(density(train$Cloud3pm,na.rm = T),col=2,main="Cloud3pm")
lines(density(train_imputed1$Cloud3pm),col=3)
plot(density(train$Cloud9am,na.rm = T),col=2,main="Cloud9am")
lines(density(train_imputed1$Cloud9am),col=3)train_imputed$Cloud3pm <- train_imputed1$Cloud3pm
train_imputed$Cloud9am <- train_imputed1$Cloud9am
test_imputed$Cloud3pm <- test_imputed1$Cloud3pm
test_imputed$Cloud9am <- test_imputed1$Cloud9amVariable Evaporation:
En esta Varible utilizaremos la media de cada estacion y con esto sustituiremos los valores faltantes.
summer <- filter(train, Season == "summer")
evaporationMeanSummer <- mean(summer$Evaporation,na.rm = TRUE)
winter <- filter(train, Season == "winter")
evaporationMeanWinter <- mean(winter$Evaporation,na.rm = TRUE)
spring <- filter(train, Season == "spring")
evaporationMeanSpring <- mean(spring$Evaporation,na.rm = TRUE)
fall <- filter(train, Season == "fall")
evaporationMeanFall <- mean(fall$Evaporation,na.rm = TRUE)
evaporationMeanSummer## [1] 7.992145
## [1] 5.715113
## [1] 4.140863
## [1] 2.447918
Aqui podemos observar como la variable de evaporacion varia bastante en funcion de la epoca del aƱo, en este caso por orden tenemos Verano, Primavera, OtoƱo e Invierno de mayor a menor, como en el resto de datos de esta variable solo utiliza un decimal haremos lo mismo para nuestras imputaciones.
evaporationMeanSummer <- round(evaporationMeanSummer,1)
evaporationMeanSpring <- round(evaporationMeanSpring,1)
evaporationMeanFall <- round(evaporationMeanFall,1)
evaporationMeanWinter <- round(evaporationMeanWinter,1)
evaporationMeanSummer## [1] 8
## [1] 5.7
## [1] 4.1
## [1] 2.4
train_imputed$Evaporation[which(is.na(train_imputed$Evaporation))] <- ifelse(train_imputed$Season == "summer", evaporationMeanSummer,
ifelse(train_imputed$Season == "winter", evaporationMeanWinter,
ifelse(train_imputed$Season == "spring", evaporationMeanSpring,
ifelse(train_imputed$Season == "fall", evaporationMeanFall,
0))))
test_imputed$Evaporation[which(is.na(test_imputed$Evaporation))] <- ifelse(test_imputed$Season == "summer", evaporationMeanSummer,
ifelse(test_imputed$Season == "winter", evaporationMeanWinter,
ifelse(test_imputed$Season == "spring", evaporationMeanSpring,
ifelse(test_imputed$Season == "fall", evaporationMeanFall,
0))))Variable Sunshine:
En esta Varible tambien utilizaremos la media de cada estacion y con esto sustituiremos los valores faltantes.
SunshineMeanSummer <- mean(summer$Sunshine,na.rm = TRUE)
SunshineMeanWinter <- mean(winter$Sunshine,na.rm = TRUE)
SunshineMeanSpring <- mean(spring$Sunshine,na.rm = TRUE)
SunshineMeanFall <- mean(fall$Sunshine,na.rm = TRUE)
SunshineMeanSummer## [1] 8.630411
## [1] 7.671015
## [1] 6.725528
## [1] 5.67726
Aqui podemos observar como la variable de Sunshine varia bastante en funcion de la epoca del aƱo, en este caso por orden tenemos Verano, Primavera, OtoƱo e Invierno de mayor a menor, como en el resto de datos de esta variable solo utiliza un decimal haremos lo mismo para nuestras imputaciones.
SunshineMeanSummer <- round(SunshineMeanSummer,1)
SunshineMeanSpring <- round(SunshineMeanSpring,1)
SunshineMeanFall <- round(SunshineMeanFall,1)
SunshineMeanWinter <- round(SunshineMeanWinter,1)
SunshineMeanSummer## [1] 8.6
## [1] 7.7
## [1] 6.7
## [1] 5.7
train_imputed$Sunshine[which(is.na(train_imputed$Sunshine))] <- ifelse(train_imputed$Season == "summer", SunshineMeanSummer,
ifelse(train_imputed$Season == "winter", SunshineMeanWinter,
ifelse(train_imputed$Season == "spring", SunshineMeanSpring,
ifelse(train_imputed$Season == "fall", SunshineMeanFall,
0))))
test_imputed$Sunshine[which(is.na(test_imputed$Sunshine))] <- ifelse(test_imputed$Season == "summer", SunshineMeanSummer,
ifelse(test_imputed$Season == "winter", SunshineMeanWinter,
ifelse(test_imputed$Season == "spring", SunshineMeanSpring,
ifelse(test_imputed$Season == "fall", SunshineMeanFall,
0))))Variables WindDir & WindSpeed: Para imputar los datos faltantes en la direccion del viento y en la velocidad utilizaremos la velocidad media que hemos obtenido antes en funcion de la direccion que tenia el aire a esa misma hora, asi nos aseguramos de dar un valor mas acorde a la realidad con los datos de los que disponemos.
Para ello primero tenemos que hacer un estudio de las variables y sus datos faltantes
trainSelect <- train %>% select(WindDir9am, WindDir3pm, WindSpeed9am,WindSpeed3pm)
gg_miss_upset(trainSelect)Con esta informacion podemos concluir que los datos que mas faltan son los de la maƱana y destacando los datos de la direccion a las 9 de la maƱana. Teniendo en cuenta que en el caso de WindDir3pm, WindSpeed9am, WindSpeed3pm, la mayoria de datos NA corresponden con valores faltantes en los demas valores no seremos capaces de imputar dichos valores con estos datos.
## # A tibble: 17 x 2
## WindDir9am avg_WindSpeed9am
## <chr> <dbl>
## 1 E 9.95
## 2 ENE 9.96
## 3 ESE 10.6
## 4 N 19.4
## 5 NE 9.89
## 6 NNE 10.5
## 7 NNW 14.6
## 8 NW 15.8
## 9 S 13.3
## 10 SE 11.2
## 11 SSE 13.4
## 12 SSW 12.7
## 13 SW 13.2
## 14 W 15.9
## 15 WNW 14.0
## 16 WSW 14.8
## 17 <NA> NA
para poder completar todos los valores vamos a imputar datos en la velocidad del viento. Para ello haremos uso del mice con el que gracias a los valores de la presion a la misma hora imputaremos los valores faltantes.
par(mfrow=c(1,1))
train_imputed0 <- mice(train_imputed[,c('Pressure9am','WindSpeed9am')], seed=2018, print = F, m = 30)
train_imputed1 <- mice::complete(train_imputed0)
test_imputed0 <- mice(test_imputed[,c('Pressure9am','WindSpeed9am')], seed=2018, print = F, m = 30)
test_imputed1 <- mice::complete(test_imputed0)
xyplot(train_imputed0, Pressure9am ~WindSpeed9am)par(mfrow=c(1,1))
train_imputed2 <- mice(train_imputed[,c('Pressure3pm','WindSpeed3pm')], seed=2018, print = F, m = 30)
train_imputed3 <- mice::complete(train_imputed2)
test_imputed2 <- mice(test_imputed[,c('Pressure3pm','WindSpeed3pm')], seed=2018, print = F, m = 30)
test_imputed3 <- mice::complete(test_imputed2)
xyplot(train_imputed2, Pressure3pm ~WindSpeed3pm)par(mfrow=c(1,2))
plot(density(train$WindSpeed9am,na.rm = T),col=2,main="WindSpeed9am")
lines(density(train_imputed1$WindSpeed9am),col=3)
plot(density(train$WindSpeed3pm,na.rm = T),col=2,main="WindSpeed3pm")
lines(density(train_imputed3$WindSpeed3pm),col=3)train_imputed$WindSpeed9am <- train_imputed1$WindSpeed9am
train_imputed$WindSpeed3pm <- train_imputed3$WindSpeed3pm
test_imputed$WindSpeed9am <- test_imputed1$WindSpeed9am
test_imputed$WindSpeed3pm <- test_imputed3$WindSpeed3pma <- data.frame(WindSpeed9am= train_imputed$WindSpeed9am, WindSpeed3pm= train_imputed$WindSpeed3pm)
train_imputed$WindGustSpeed[which(is.na(train_imputed$WindGustSpeed))] <- rowMeans(a, na.rm=TRUE)
b <- data.frame(WindSpeed9am= test_imputed$WindSpeed9am, WindSpeed3pm= test_imputed$WindSpeed3pm)
test_imputed$WindGustSpeed[which(is.na(train_imputed$WindGustSpeed))] <- rowMeans(b, na.rm=TRUE)Transformaciones de variables cualitativas
train_imputed[c("Temp9am_imp","Temp3pm_imp" ,"Pressure9am_imp","Pressure3pm_imp","MaxTemp_imp","MinTemp_imp", "Humidity9am_imp","Humidity3pm_imp","Rainfall_imp")] <- list(NULL)
test_imputed[c("Temp9am_imp","Temp3pm_imp" ,"Pressure9am_imp","Pressure3pm_imp","MaxTemp_imp","MinTemp_imp", "Humidity9am_imp","Humidity3pm_imp","Rainfall_imp")] <- list(NULL)train_imputed <- dummy_cols(train_imputed, select_columns = c("Location", "Season","WindGustDir","WindDir9am","WindDir3pm"))
test_imputed <- dummy_cols(test_imputed, select_columns = c("Location", "Season","WindGustDir","WindDir9am","WindDir3pm"))Transformaciones de variables cuantitativas La única que no parece seguir una distribución normal es la variable RainFall
p1 <- train_imputed %>% select(Rainfall, Season) %>%
na.omit() %>%
ggplot(aes(x=Rainfall, colour=Season)) +
geom_density()
#(train_imputed$Rainfall)+0.01
aux <- filter(train_imputed,train_imputed$Rainfall>0)
p2 <- aux %>% mutate(inv_Rainfall = 1/ (Rainfall)) %>%
select(Season, inv_Rainfall) %>%
na.omit() %>%
ggplot(aes(x=inv_Rainfall, colour=Season)) +
geom_density()
p3 <- train_imputed %>% mutate(log10_Rainfall = log10(Rainfall + 0.01)) %>% #+ 1.01
select(Season, log10_Rainfall) %>%
na.omit() %>%
ggplot(aes(x=log10_Rainfall, colour=Season)) +
geom_density()
hist(train_imputed$Rainfall)Con esta informacion podemos concluir que para transformarla lo mejor que podemos hacer es utilizar la inversa sin tener en cuenta los ceros de la variable rainfall, que equivaldrian a un āNoā en la variable RainToday
#Resto de variables: siguen distribución normal
qqnorm(train_imputed$MinTemp, ylab="MinTemp")
qqline(train_imputed$MinTemp, col="red")## [1] 9824 88
## [1] 1228 88
## [,1]
## Date 0.00000000
## Season 0.00000000
## Location 0.00000000
## MinTemp 0.00000000
## MaxTemp 0.00000000
## Rainfall 0.00000000
## Evaporation 0.00000000
## Sunshine 0.00000000
## WindGustDir 11.24796417
## WindGustSpeed 0.00000000
## WindDir9am 7.21701954
## WindDir3pm 2.33102606
## WindSpeed9am 0.00000000
## WindSpeed3pm 0.00000000
## Humidity9am 0.00000000
## Humidity3pm 0.00000000
## Pressure9am 0.00000000
## Pressure3pm 0.00000000
## Cloud9am 0.00000000
## Cloud3pm 0.00000000
## Temp9am 0.24429967
## Temp3pm 0.17304560
## RainToday 0.00000000
## RISK_MM 0.00000000
## RainTomorrow 0.00000000
## AvgTemp 0.04071661
## Temp 0.02035831
## Pressure 1.81188925
## Humidity 0.09161238
## Location_Sydney 0.00000000
## Location_Adelaide 0.00000000
## Location_Melbourne 0.00000000
## Location_Canberra 0.00000000
## Season_summer 0.00000000
## Season_fall 0.00000000
## Season_spring 0.00000000
## Season_winter 0.00000000
## WindGustDir_SSE 11.24796417
## WindGustDir_ESE 11.24796417
## WindGustDir_S 11.24796417
## WindGustDir_SW 11.24796417
## WindGustDir_WSW 11.24796417
## WindGustDir_NW 11.24796417
## WindGustDir_N 11.24796417
## WindGustDir_NNE 11.24796417
## WindGustDir_NNW 11.24796417
## WindGustDir_WNW 11.24796417
## WindGustDir_W 11.24796417
## WindGustDir_NA 0.00000000
## WindGustDir_NE 11.24796417
## WindGustDir_SE 11.24796417
## WindGustDir_SSW 11.24796417
## WindGustDir_ENE 11.24796417
## WindGustDir_E 11.24796417
## WindDir9am_SSW 7.21701954
## WindDir9am_E 7.21701954
## WindDir9am_SE 7.21701954
## WindDir9am_WNW 7.21701954
## WindDir9am_W 7.21701954
## WindDir9am_NE 7.21701954
## WindDir9am_NW 7.21701954
## WindDir9am_SW 7.21701954
## WindDir9am_NA 0.00000000
## WindDir9am_N 7.21701954
## WindDir9am_SSE 7.21701954
## WindDir9am_S 7.21701954
## WindDir9am_ENE 7.21701954
## WindDir9am_WSW 7.21701954
## WindDir9am_NNW 7.21701954
## WindDir9am_NNE 7.21701954
## WindDir9am_ESE 7.21701954
## WindDir3pm_S 2.33102606
## WindDir3pm_E 2.33102606
## WindDir3pm_SW 2.33102606
## WindDir3pm_W 2.33102606
## WindDir3pm_ESE 2.33102606
## WindDir3pm_WNW 2.33102606
## WindDir3pm_N 2.33102606
## WindDir3pm_SSW 2.33102606
## WindDir3pm_SE 2.33102606
## WindDir3pm_NNW 2.33102606
## WindDir3pm_ENE 2.33102606
## WindDir3pm_NW 2.33102606
## WindDir3pm_WSW 2.33102606
## WindDir3pm_NA 0.00000000
## WindDir3pm_NE 2.33102606
## WindDir3pm_NNE 2.33102606
## WindDir3pm_SSE 2.33102606
El modelo utilizado es una regresion logistica, sacamos los valores de la recta para cada variable mediante lasso
Aplicamos la estandarización con las medias y varianzas de train.
#Estandarización
numeric_vars <- c("MinTemp","MaxTemp", "Rainfall", "Evaporation", "Sunshine", "WindGustSpeed",
"WindSpeed9am", "WindSpeed3pm", "Humidity9am", "Humidity3pm","Pressure9am" ,
"Pressure3pm" , "Temp9am", "Temp3pm", "RISK_MM")
preprocessParams <- preProcess(train_imputed[numeric_vars], method=c("center", "scale"))
train_imputed[numeric_vars] <- predict(preprocessParams, train_imputed[numeric_vars])
#Aplicamos la estandarización con las medias y varianzas de train
test_imputed[numeric_vars] <- predict(preprocessParams, test_imputed[numeric_vars])Creación variables dummies para Location, Season, WindGustDir y borrado de las anteriores.
#Procesado de variables cualitativas: creación variables dummies
train_imputed <- dummy_cols(train_imputed, select_columns = c("Location", 'Season', 'WindGustDir'))
test_imputed <- dummy_cols(test_imputed, select_columns = c("Location", 'Season', 'WindGustDir'))
train_imputed$Location <- NULL
train_imputed$Season <- NULL
train_imputed$Date <- NULL
train_imputed$WindGustDir <- NULL
test_imputed$Location <- NULL
test_imputed$Season <- NULL
test_imputed$Date <- NULL
test_imputed$WindGustDir <- NULL
test_imputed$WindDir3pm_NA <- NULL
test_imputed$WindDir9am_NA <- NULLtrain_imputed_sinNA = na.omit(train_imputed)
x_train = train_imputed_sinNA[, !names(train_imputed_sinNA) %in% c("RainTomorrow")]
y_train = train_imputed_sinNA$RainTomorrow
test_imputed_sinNA = na.omit(test_imputed)
x_test = test_imputed_sinNA[, !names(test_imputed_sinNA) %in% c("RainTomorrow")]
y_test = test_imputed_sinNA$RainTomorrow
dim(test_imputed_sinNA)## [1] 1016 82
## [1] 8169 84
#Primero vemos Lasso
x = model.matrix(RainTomorrow~ ., train_imputed_sinNA)[,-1]
lambda_seq <- 10^seq(2, -2, by = -.1)
cv.out <- cv.glmnet(x, y_train, alpha = 1, lambda = lambda_seq)
plot(cv.out)# identifying best lamda
best_lam <- cv.out$lambda.min
lasso_best <- glmnet(x, y_train, alpha = 1, lambda = best_lam)
c<-coef(lasso_best,s=best_lam,exact=TRUE)
inds<-which(c!=0)
variables<-row.names(c)[inds]
variables## [1] "(Intercept)" "Sunshine" "WindGustSpeed"
## [4] "WindDir9amN" "WindDir9amS" "WindDir9amSSW"
## [7] "WindDir3pmNNE" "WindSpeed3pm" "Humidity3pm"
## [10] "Pressure3pm" "Cloud3pm" "RainToday"
## [13] "RISK_MM" "Location_Sydney" "Location_Adelaide"
## [16] "Season_summer" "Season_winter" "WindGustDir_NNW"
## [19] "WindGustDir_ENE" "WindDir9am_SSW" "WindDir9am_S"
## [22] "WindDir9am_NNE" "WindDir3pm_NNE"
new_train = train_imputed_sinNA %>% select("Sunshine","WindGustSpeed", "WindDir9am_N", "WindDir9am_NNE", "WindDir9am_S",
"WindDir3pm_NNE", "Humidity3pm", "Pressure3pm", "Cloud3pm", "RainToday",
"Location_Adelaide", "Location_Sydney", "Season_summer","WindGustDir_ENE", "WindGustDir_NNW",
"WindDir9am_N", "WindDir9am_NNE", "WindDir9am_S", "WindDir3pm_NNE", "RainTomorrow")
new_test = test_imputed_sinNA %>% select("Sunshine","WindGustSpeed", "WindDir9am_N", "WindDir9am_NNE", "WindDir9am_S",
"WindDir3pm_NNE", "Humidity3pm", "Pressure3pm", "Cloud3pm", "RainToday",
"Location_Adelaide", "Location_Sydney", "Season_summer","WindGustDir_ENE", "WindGustDir_NNW",
"WindDir9am_N", "WindDir9am_NNE", "WindDir9am_S", "WindDir3pm_NNE", "RainTomorrow")
#Modelo
#Entrenamiento
glm_model_train = glm(RainTomorrow~ ., data=new_train, family= binomial)
#Test
glm_test = predict(glm_model_train, newdata = new_test, type = "response")
#Evaluación modelo
#Summary
summary(glm_model_train)##
## Call:
## glm(formula = RainTomorrow ~ ., family = binomial, data = new_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6951 -0.5643 -0.3426 -0.1585 2.8729
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.39082 0.10181 -23.482 < 2e-16 ***
## Sunshine -0.50184 0.03863 -12.990 < 2e-16 ***
## WindGustSpeed 0.47727 0.04124 11.573 < 2e-16 ***
## WindDir9am_N 0.09915 0.09587 1.034 0.301033
## WindDir9am_NNE 0.39544 0.13169 3.003 0.002675 **
## WindDir9am_S -0.33361 0.14500 -2.301 0.021408 *
## WindDir3pm_NNE 0.49836 0.17700 2.816 0.004867 **
## Humidity3pm 0.91946 0.04239 21.690 < 2e-16 ***
## Pressure3pm -0.51777 0.03936 -13.153 < 2e-16 ***
## Cloud3pm 0.05060 0.01401 3.611 0.000306 ***
## RainToday 0.55232 0.07221 7.648 2.03e-14 ***
## Location_Adelaide 0.63441 0.08212 7.725 1.12e-14 ***
## Location_Sydney 0.22862 0.08438 2.709 0.006743 **
## Season_summer -0.16318 0.08425 -1.937 0.052763 .
## WindGustDir_ENE -0.45050 0.18304 -2.461 0.013848 *
## WindGustDir_NNW 0.34857 0.12695 2.746 0.006040 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8823.1 on 8168 degrees of freedom
## Residual deviance: 6117.7 on 8153 degrees of freedom
## AIC: 6149.7
##
## Number of Fisher Scoring iterations: 5
#Tabla de ganancia
logistic_gains_table <- blr_gains_table(glm_model_train, data = new_train)
#Curva ROC
blr_roc_curve(logistic_gains_table)#Matriz de confusión
umbral_dec = 0.46
glm_test <- ifelse(glm_test >= umbral_dec, 1, 0)
glm_test <- factor(glm_test, levels = c(0,1))
tabla_conf <- table(glm_test, new_test$RainTomorrow)
caret::confusionMatrix(tabla_conf, positive = '1')## Confusion Matrix and Statistics
##
##
## glm_test 0 1
## 0 731 115
## 1 59 111
##
## Accuracy : 0.8287
## 95% CI : (0.8041, 0.8514)
## No Information Rate : 0.7776
## P-Value [Acc > NIR] : 3.216e-05
##
## Kappa : 0.4569
##
## Mcnemar's Test P-Value : 3.052e-05
##
## Sensitivity : 0.4912
## Specificity : 0.9253
## Pos Pred Value : 0.6529
## Neg Pred Value : 0.8641
## Prevalence : 0.2224
## Detection Rate : 0.1093
## Detection Prevalence : 0.1673
## Balanced Accuracy : 0.7082
##
## 'Positive' Class : 1
##